library(MASS)
library(anchors)

setwd("C:/XunCao/Research/Space_RA/Maps/Bosnia/bosnia survey and gis/survey")
dat<-dget("wobc_052007")
dat<-cbind(dat, rep(1, dim(dat)[1]))
colnames(dat)<-c(colnames(dat)[-dim(dat)[2]], "ones")

dat.no.na<-subset(dat, dat$nat.trust!=9
                                 &dat$vladimir!=9
                                 &dat$ethnicrelations!=9
                                 &dat$currentsit!=9
                                 &dat$ethnicfriends!=9
                                 &dat$closestfriends!=9
                                 &dat$pride!=9
                                 &dat$organization!=9
                                 &dat$violence!=9
                                 &dat$materialstatus!=9
                                 &dat$education!=9)


dim(dat.no.na)

defdta <- list(z1 = c("vladimir"), 
               y0 = c("nat.trust"),                      
               x0 = c("age", "gender", "education", "materialstatus","pride", 
                      "currentsit", "ethnicfriends", "closestfriends", "ethnicrelations", "violence", "organization"
                      ,
                      "avar" , "darghinn", "kabardinian", "ossetin",  "russian",  
                      "bosniac", "serb", "croat"),     # covariates for mu      
               v0 = c("age", "gender", "education", "materialstatus","pride", 
                      "currentsit", "ethnicfriends", "closestfriends", "ethnicrelations", "violence", "organization"
                      ,
                      "avar" , "darghinn", "kabardinian", "ossetin",  "russian",  
                      "bosniac", "serb", "croat",
                      "ones")                          # covariates for gammas
               )
out.chopit  <- chopit(dat.no.na, defdta, use.linear=T) # use 

summary.chopit(out.chopit) 



####

n.sim <- 100
n.parm<-out.chopit$chopit.parm$nvec

hess<-solve(out.chopit$chopit.hess) # this is the covariance matrix
colnames(hess)<-rownames(hess)<-n.parm
theta.gamma <- out.chopit$chopit.parm$pvec # mean of all coefficients
theta.gamma.h <- hess
dimnames(theta.gamma.h)

set.seed(391)
sim <- mvrnorm(n=n.sim, mu=theta.gamma, Sigma=theta.gamma.h, tol=1e-6)


n.cat  <- 5       ## number of choice categories per question
n.vign <- 1       ## number of vignettes quesitons
n.self <- 1       ## number of self questions



# now pick different nationalities:

pick<-c("Avar", "Dargin", "Kabardin", "Ossetian", "Russian",  "Bosniac", "Serb", "Croat")
nat<-c("avar" , "darghinn", "kabardinian", "ossetin", "russian", "bosniac", "serb", "croat")
dta  <- dat.no.na

ylim <- c(0,9)
xlim <- c(-3.5,5)

plot(xlim,ylim,type="n",xlab="",ylab="",axes=FALSE)#, main=mn)National Trust: Latent Variable Values
axis(1);  axis(2, labels=pick, at=c(1:8), las=1, padj=0, hadj=0, col="lightgray", tick=F)


## have to do for each nationality: I know it is a pain ... but 

h<-1 # pick the group, say if it is 1, then it is avar ... 
bdta  <- dta[ dta$avar == 1, ]  ### make dta$avar == 1 ==> only choose Avar 

n <- NROW(dta)

bx0 <- as.matrix((bdta[,defdta$x0]))# b for Bosnia
bv0 <- as.matrix((bdta[,defdta$v0]))

txt0 <- dimnames(theta.gamma.h)[[1]][(length(defdta$v0)*4+3):(length(defdta$v0)*4+length(defdta$x0)+2)]
txt1 <- dimnames(theta.gamma.h)[[1]][length(defdta$v0)*4+2]
txt2 <- dimnames(theta.gamma.h)[[1]][1:(length(defdta$v0)*4)]
txt3 <- dimnames(theta.gamma.h)[[1]][length(defdta$v0)*4+1]

se.vign1 <- exp(as.real(sim[,txt3]))

bmu <- nmu <- jtheta <- btaus <- ntaus <- NULL
for (i in 1:n.sim) {
 se <- se.vign1[i]
 sim.beta  <- sim[i,txt0]
 sim.theta <- sim[i,txt1]
 sim.gamma <- matrix(sim[i,txt2],nrow=NCOL(bv0))

 bmu  <- rbind(bmu, bx0 %*% sim.beta)
 bvg  <- bv0 %*% sim.gamma ##
 btau <- bvg # exp(bvg)          ## exponentiate all
 btau[,1] <- bvg[,1]       ## re-make first tau into  levels
 #taus <- row.cumsum(btau) ## get cumulative
 taus<-btau
 taus[,2]<-taus[,1]+taus[,2]
 taus[,3]<-taus[,2]+taus[,3]
 taus[,4]<-taus[,3]+taus[,4]
 btaus <- rbind(btaus,taus)

 jtheta<- rbind(jtheta,sim.theta)
}

# plot the latent variable for the DV: 95% CI
x<-c(mean(bmu)-1.96, mean(bmu)-1.96, mean(bmu)+1.96, mean(bmu)+1.96)
y<-c(h+.2, h-.2, h-.2, h+.2)
polygon( x, y, col="grey",border="grey")

# plot the distribution of thresholds:
hl<-c(-.2, .10, -.1, .2)
for (i in 1:(n.cat-1)) {sum.post<-c(quantile(btaus[,i], c(.025, .975))[1], mean(btaus[,i]), quantile(btaus[,i], c(.025, .975))[2])
                        segments(sum.post[1], h+hl[i], sum.post[3], h+hl[i], col="black", lwd=3); points(sum.post[2], h+hl[i], pch=19, cex=1.5)}



### finally, the ropeladder for 95% CI of coefs:

va<-c("Age", "Gender", "Education", "Material status", "Pride", "Current situation",
      "Ethnic friends", "Closest friends", "Ethnic relations", "Violence", "Organization", 
      "Avar", "Dargin", "Kabardin", "Ossetian", "Russian", "Bosniak", "Serb", "Croat")


op.m<-c(0.001, 0.000, 0.088, 0.020, 0.044, -0.059, -0.140, 0.253, -0.159, 0.151, 0.122, 0.070, 0.020, 0.519, 0.032, -0.154, -0.694, -.611, -0.690)
op.se<-c(0.001, 0.038, 0.018, 0.028, 0.022, 0.036,  0.023,  0.027,  0.035, 0.045, 0.103, 0.089, 0.108, 0.131, 0.121,  0.062,  0.069, 0.068, 0.083)


chop.m<-c(-0.001, -0.063, 0.094, -0.016, 0.112, -0.117, -0.022, 0.308, -0.062, 0.251, -0.083, -0.001, 0.438, 0.670, 0.419, 0.045, -0.734, -0.648, -0.382)
chop.se<-c(0.002, 0.062, 0.029, 0.044, 0.034, 0.058, 0.033, 0.043, 0.057, 0.073, 0.166, 0.159, 0.176, 0.175, 0.183, 0.103, 0.114, 0.117, 0.137)


par(mfrow=c(1,2))
ylim <- c(1,19)
xlim <- c(-1.2,1)

plot(xlim,ylim,type="n",xlab="",ylab="",axes=FALSE, main="ordered probit");abline(v=0, col="gray")
axis(1);  axis(2, labels=va, at=c(c(19:1)-.1), las=1, padj=0, hadj=0, col="lightgray", tick=F)
segments(op.m-1.96*op.se, 19:1, op.m+1.96*op.se, 19:1, col="gray", lwd=2); points(op.m, 19:1, pch=19, cex=1)

plot(xlim,ylim,type="n",xlab="",ylab="",axes=FALSE, main="chopit");abline(v=0, col="gray")
axis(1);   axis(2, labels=va, at=c(c(19:1)-.1), las=1, padj=0, hadj=0, col="lightgray", tick=F)
segments(chop.m-1.96*chop.se, 19:1, chop.m+1.96*chop.se, 19:1, col="gray", lwd=2); points(chop.m, 19:1, pch=19, cex=1)
